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ABSTRACT 

This thesis studies the roll loop of an operational 
surface to air missile as an example of a typical self 
adaptive control system. The missile studied also exhibits 
a phenomenum called ROLL WANDER, a low frequency oscillation 
in missile roll orientation. Because of this, the thesis is 
developed in three major areas of study. First, is the study 
of the missile roll control loop. Second, is a study of the 
phenomenon of ROLL WANDER, with the production of this pheno- 
menon in the computer simulation. Third and finally is the 
implementation, in the computer simulation, of a self adaptive 
notch filter to suppress the occurence of ROLL WANDER. 
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I. CONTROL SYSTEMS 



A . INTRODUCTION 

The objective of this thesis was to study a self- 
adaptive control system through the use of a computer 
simulation. The control system studied was the roll control 
loop of an operational surface to air missile. 

This particular control loop uses a bistable actuator 
(a nonlinear element) as the primary mover of the loop. 

This method of control, a full throw "bang bang" system is 
commonly used in control systems. The control loop is driven 
either full positive or full negative, depending upon the 
polarity of the error signal. Since the system reacts with 
full response, the output is driven rapidly in opposition to 
the error and the error signal is eventually driven to the 
opposite polarity. The control then changes sign and the 
output is driven in the opposite direction. In this manner, 
the error signal is nulled out. 

Bang bang control systems are in wide-spread use for such 
systems as surface to surface missiles, surface to air 
missiles, and some types of surface and submarine torpedoes. 
All of these systems require an accurate, reliable type one 
system producable at a minimum of expense. A full throw 
control system can provide these requirements. 

In the case of a surface to air or air to air missile, 
however, there are added problems that the control system 
must be able to cope with. It must be able to operate at 



9 



various altitudes and at various speeds. The correct control 
response to an error signal of a given magnitude at high 
altitude (thin atmosphere) would be a drastic over response 
to the same error signal at low altitude (dense atmosphere). 

To correct for this, the control system must be self adaptive. 
It must respond to the error signal at a rate appropriate 
with its medium. Another constraint placed upon the control 
system is the requirement not to disturb the missile guidance 
and seeking/tracking mechanisms. While a full throw control 
system is inexpensive and reliable, it does introduce a 
constant oscillation in the position output. This can be 
disruptive to the other missile systems mentioned. 

The surface to air guided missile studied in this thesis 
utilizes a bang bang type of control system which is both 
self-adaptive and non-disruptive to other missile functions. 

B. DITHERING IN THE ROLL SYSTEM 

Figure 1-1 is a block diagram of the roll control loop. 

The bistable element provides the control command to the 
tail control surface actuators. The use of this nonlinear 
element in the control loop causes a limit cycle oscillation 
to appear. It is this limit cycle oscillation which provides 
the self adaptive feature of the autopilot, and at the same 
time allows for smooth loop response while using a bang bang 
type control system. With no error or command signal present, 
the limit cycle oscillation provides a constant amplitude 
dither signal to the tail control surface servos. The dither 
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with zero error signal has equal negative and positive 
components; and as a result, the clockwise and counter 
clockwise roll torques produced by the tail surfaces cancel. 

The missile therefore dithers slightly in roll but maintains 
a constant average roll orientation. The amplitude of the 
dither signal, and thus the adaptive gain characteristics of 
the autopilot is a function of G, the roll moment coeffi- 
cient. G is a measure of the aerodynamic gain, a property 
of the medium in which the missile is traveling 

C. SELF-ADAPTIVE FEATURE 

The autopilot derives its self-adaptive feature by 
varying the amplitude of limit cycle oscillation, or dither, 
in the roll control loop. A study of Figure 1-1 shows that 
there are three basic loops. The outer two loops provide 
position (<f>') and velocity (<J>') feedback paths for stability 
and damping. The inner loop also provides some rate (<J>’) 
feedback, but of greater importance, the inner loop coupled 
with the bistable, element produces the limit cycle oscillation. 
The roll loop of the missile was used for the dithering loop, 
since dithering in the other control loops would not be 
practical. Dithering in the pitch loop would be intolerable 
due to the amount of body motion sensed by the head rate 
gyros. Similarly, a dither in the squeeze mode would not 
provide the aerodynamic adaptive gain features. 

Equation 1 is the dynamic equation for the rigid body 
motion of the missile. 
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( 1 ) 



(j) = g s - r (J> 



o°. <j) = -S_ (la) 

S (s-F) 



where : (j- 


= the aerodynamic moment 


V 


= the aerodynamic damping 



This transfer function, called the "aerodynamic roll 
transfer function" produces the inherent self-adaptive 



feature . 




G = 


WBIASq' n 2 ^3C LS (lb) 

-S-R 


where A 


= static air pressure ratio 


n 


= mach no. 



Ik = roll inertia 

roll damping coefficient 



F= 


1431 ASdn^TAbC^ do 


where: \J 


= velocity 


S,d,b 


= aerodynamic constants 



As can be seen by Equations (lb) and (lc), the gain of 
the roll transfer function is dependent on flight conditions. 
Therefore, as the missile’s environment changes, so does the 
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gain of the roll transfer function and, accordingly, the 
amplitude of the position loop's limit cycle oscillation, or 
dither. 

D. ROLL CONTROL 

The controlling action of the self-adaptive loop is 
shown in Figure 1-2. The tail control surface drive motors 
are driven by the output of the bistable element at a constant 
rate in either the positive or negative direction. When an 
error signal is present, it appears as a dc bias signal 
summed with the roll loop's natural dither oscillation 
(Figure l-2a). The algebraic sum of the two signals is fed 
to the bistable element (Figure l-2b). As its name implies, 
the bistable element has only two states, full positive and 
full negative. When the value of K^F^ (the summation signal) 
passes through zero the bistable element shifts states. The 
action of the dc bias produced by the error signal is to 
cause the square wave output of the bistable element to 
remain positive or negative longer during a complete cycle, 
as appropriate (Figure (l-2c). Since the tail control 
surface motors act to integrate this output, the result is 
a dithering movement of the control surfaces modulated by 
the appropriate control surface movement to null the roll 
error signal. 

As previously discussed, the gain of the roll transfer 
function determines the amplitude of the dither oscillation. 

As a missile flies higher, the gain decreases and therefore 
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the amplitude of the dither decreases. The effect of the 
dc bias is more pronounced with the result that the control 
loop is more sensitive and has higher gain. Conversely, 
at lower altitudes or higher speeds, the dither amplitude 
is greater and the error signal has less effect. The loop 
gain is reduced (Figure l-3a) . 

The dither imposes a maximum signal strength allowable 
for the error signal. If at a given dither amplitude an 
error signal is encountered which is greater than the dither 
amplitude, then the oscillation is shifted completely to 
one side- of the axis (Figure l-3b). The bistable element 
stays continuously in one state and the tail control surfaces 
are driven continually in one direction. Loop response 
diverges and the loop becomes unstable (Figure l-3c) . To 
prevent this occurrence, the error signal is electronically 
limited. 

Smooth response in the roll loop is obtained by super- 
imposing the dc error signal on the dither. In fact the 
error signal is not a pure dc signal, but is also a time 
varying waveform. It is important to remember, therefore, 
that the dither frequency must be sufficiently higher than 
the missile control frequencies (error signal frequencies) 
so that the gain adaptive feature can be maintained. This 
imposes a severe lower bound on an acceptable dither 
frequency . 
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E. DESCRIPTION OF LOOP COMPONENTS 



Figure 1-1 can be divided into three major portions: 
the tail control surface actuators, the air frame, and the 
autopilot electronics package with associated sensors. 

The original design work on the missile started with a 
given set of actuators and the air frame's equations of 
motion. The autopilot control equations were then written 
to drive the air frame with the given hardware . 

Three major parameters are input or sensed within the 
roll control system. They are <j> c , the commanded roll angle, 

K^Fg, the bistable element forcing function (commonly called 

« 

"delta dot called for") and r , the tail roll rate as 

« 

measured by the roll rate gyro. <J>' is integrated by the roll 
free gyro to produce the sensed missile roll angle <J>*. 4> 

and <J>' are then differenced to obtain <J> , the roll error 
signal . 

1 . Bistable Forcing Function 



Equation -2 is the equation of the basic roll forcing 



function into the bistable switcher. 



~ L'fb + 

*K<V 




where : Fs = roll control signal 




\ = the electronically 

• T f* i - o onnnn c? rr-n Q 



the electronically limited error signal. 
If the error signal is allowed to get 
too large, the loop will stop dithering 
and diverge 
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(|/ = roll rate output of the roll rate gyro 

K$,K$ are system constants 

Preliminary analysis showed that high frequency 
resonances caused by the missile's elastic body bending mode 
were possible. To suppress signals at resonant frequencies 
and at the same time maintain the dither frequency, a 

cascaded set of high frequency (w^-Q^er > w dither^ notch 
filters were added. The final form of the roll forcing 
function (Equation 2) becomes: 




where 




(3a) 



(3b) 



(3c) 
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2 . 



Roll Rate Gyro 



The transfer function for the roll rate gyro is 




where : 



- / 

<p = the gyro output signal 

<fi = missile roll rate 

f = mechanical damping coefficient 



( 4 ) 



3 . Roll Free Gyro 

The roll free gyro has a transfer function of unity 
and is considered linear over all normal operating ranges . 

^ . Tail Servo Transfer Function 

The tail servo transfer function (using a rigid body 
model vice an elastic body model) can be described as 



a _ 



& 



e 



(t-+ at)s 

s 



( 5 ) 



where: T = the actuator time delay 

AT 

= an approximation to the extra phase lag intro- 
duced by high frequency noise filtering, switch 
hystersis and aerodynamic plane lag due to the 
tail dithering. 

5 . Tail Surface Hinge Moment 



Shown as a positive feedback path around the tail 
servo transfer function is K^K^. represents the hinge 

moment due to tail position. Aerodynamic loading forces 
acting on the tail control surfaces produce this hinge move- 
ment torque at the output of the tail actuators. The loop 
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is shown as a positive feedback loop, since the greater the 
effective tail angle, the greater the hinge moment torque 
produced; which, in turn, tends to increase the effective 
tail angle. is a conversion factor, and is given as a 

system constant. 

In this thesis, the effect of hinge moment relative 
to other system parameters was negligible; and in all 
simulations it was set to zero. Hinge moment was included 
in the model for completeness. 

6 . F and F 

a squeeze 



F is the pitch control signal. It is produced in 

cl 

the pitch control loop, then summed with the roll control 
signal as an input into the bistable element. F is 

the control signal generated by the autopilot to compensate 
for the misalignment of the tail surfaces. Misalignment, 
resulting from mechanical tolerances and system slop in each 
of the four tail surfaces causes torques on the missile 
frame; counteracting torques may not produce roll or pitch, 
but will produce frame stress or squeeze. 

7 . Biases 

The effect of biases in the roll loop is to cause 
the missile to assume a roll attitude other than that estab- 
lished at time of launch. Maintenance of roll orientation 
is essential to the guidance portion of the missile. The 
various components of the roll loops in addition to the pitch 
loop through F and the squeeze signal, all 

introduce biases. 
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These loop biases are counteracted by the filter 

block: 

Rs)=(-^-) (6 

The equivalent roll error of each contributor can be deter- 
mined by solving for <p in the following simplified static 
roll control equation (7), with each bias substituted in 
terms of deg or deg/sec in the appropriate term. 



si= k u 




.33lKo<P 



m; 



bw5 



( 7 ) 
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II. MODELING THE MISSILE 



A. INTRODUCTION 

There were two major areas of study in the research. 

First was the construction and testing of an accurate 
computer model to simulate the missile. Second was an 
investigation of an abnormal phenomenon experienced in the 
missile roll control loop called "roll wander." An explana- 
tion of this problem and a proposed solution will be 
presented in subsequent chapters. 

B. SIMULATION 

Due to the complexities of the missile control system 
and the inherent nonlinear operation introduced by the 
bistable element, it was decided to model the control system 
using a completely digital simulation. The model was there- 
fore constructed using IBM’s Digital Simulation Language 
(DSL). DSL was chosen for the model because under its format 
complex transfer functions could be described by single 
blocks with accompanying data statements. This eliminated 
the need for writing state variable equations for as high as 
eighth order systems, and thereby significantly reduced the 
amount of coding required to model the missile. A listing 
of the DSL program with a more detailed commentary on the 
operation of the DSL program is provided in Appendices A 
and B. 
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Shown in Figure 2-1 is the missile roll control loop as 
it was arranged for encoding into DSL. The letter names 
for the different function blocks correspond to the variable 
names used in the simulation program. The DSL function TRNFR 
used to represent the various transfer functions imposes one 
restriction. It requires that the polynomial in the denomi- 
nator be of a greater order than the polynomial in the 
numerator. The reason for this is because the TRNFR block 
determines the correct output response by solving the 
corresponding system of first order differential equations. 
This requires the polynomial in the denominator to have the 
highest order to allow for the construction of the correct 
first order system. To accommodate this restriction, certain 
blocks (B and HP) were formed by combining other blocks (as 
seen by a comparison of Figure 1-1 and Figure 2-1) . This 
altered the appearance of the block diagram, but not the 
behavior of the loop. 

Block D represents the mechanical system’s response time 
to the bistable element’s output signal. A finite time 
occurs between the electrical command and the mechanical 
movement of the tail control surfaces . This was modeled in 
the computer program by the DSL function block DELAY. DELAY 
is a dead-time function. Sample points of the input expression 
are saved at specified intervals. The output is set equal to 
zero until the delay interval expires at which time the output 
has the value of the input at the start of the delay. If 
necessary, linear interpolation is used between sample points. 
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SIMULATION LOOP 
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This ability, gained using an all digital simulation to 
accurately and easily adjust the delay interval, was one of 
the major factors in choosing the digital simulation over an 
equivalent analog simulation. 

The bistable element was modeled using a standard 
FORTRAN logical IF statement. The bistable element has no 
zero point, only full positive and full negative. The 
logical IF was written so that when K^F^ was greater than 
zero, the bistable output was positive. When K^F^ was zero 
or negative, the bistable output was negative. 

The aerodynamic gain constants G and F were calculated 
at the beginning of each run. The parameters necessary to 
calculate these values and the parameters specifying the 
other systems constants were specified at the beginning of 
a given run. These constants could be changed during or at 
the end of the run, thus giving the capability of multiple 
runs within one DSL job. 

Difficulty was initially encountered due to the step size 
used. The TRNFR blocks are comprised of a series of inte- 
grators simulated by DSL by a call to an integration routine 
of a type specified by the user. To conserve computer time, 
a fixed interval fourth order Runga Kutta method was used. 
However, it is a property of this integration routine as 
well as many other numerical integration routines, that if 
the frequency components of the function to be integrated 
approach the frequency of the integration routine, specified 
by the routine step size, the integration method produces a 
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divergent error and becomes unstable. Such was the case with 
this program. Initially, a step size of 0.005 was used. 
Problems were encountered with the inclusion of the high 
frequency notch filters. The notch filter represented by a 
denominator of: 



has a center frequency of 108 hz. Using a step size with 
an integration frequency of only 200 hz was insufficient, 
and caused divergence. 

Optimal step size was obtained by testing the loop using 
a variable step size Runga Kutta method with different error 
limits specified. By specifying a tolerance of 0.001 
absolute error an average integration step size of 0.0004 
was obtained. Thus, a frequency of 2500 hz was sufficiently 
high to avoid numerical integration problems. During the 
actual runs, a final step size of 0.0001 was used. 

C. TESTING THE MODEL 

Except for the problem of determining the program inte- 
gration step size, the only difficulties encountered with 




( 8 ) 



has a center frequency of 159 hz. 



Similarly, the notch filter represented by: 




(9) 
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the initial simulation runs were those of writing the standard 
program logic and debugging routine syntactical errors. Once 
constructed, the model was tested under normal operating 
conditions to ensure that it behaved the same as the real 
missile. The model was allowed to "fly" with no input 
(Figure 2-2). The initial jump noted is the transient 
behavior of the loop. Since the parameters of the loop were 
continually being changed from one run to the next, initial 
conditions of the integrators were not calculated; rather all 
integrators were set initially to zero, and the transient 
response accepted. The low amplitude oscillation is the 
bistable dithering action. The dither appeared in the 
position loop because of the step size of integration used. 

By decreasing step size, the amplitude of the dither could 
be reduced. However, decreasing step size proportionately 
increased run time. A compromise was reached using a step 
size of 0.0001. 

The closed loop transfer function for a roll command 
( 4>/4> ) obtained by writing the loop equation for Figure 1-1, 
(with K = 0) is: 




Kq> S 3 





*^s(i+ r f b ) 

It- 

1 ^ eooc 



Gg_ Gr. 



4 > 



4 - 




( 10 ) 
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where: = transfer function of the roll rate gyro 

- torsional notch filter 

Applying the final value theorem for a step input 
(4> c = 4^/s), to Equation (10) gives: 

L'.m <j) _ <j) t K t K t e 1 -*■•%) 

Similarly, the closed loop expression for a tail induced 
torque (4> c /6 c ) , also obtained from Figure 1-1 is: 

QO O 

Sc “Sc 




(ID 



( 12 ) 
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Then for a step input (6 c = const/s), and applying the final 
value theorem to Equation (12) 

• • 

Lim — O (13) 

"( — i- 00 

The model was thus tested using a step input of -45° for 

4> and a step input of 3° for 6 . The response to the step 
c c 

input of <p is shown in Figure (2-3) • For a roll command or 

a roll error signal, the missile will roll and maintain the 

desired angle. The response of the model compared favorably 

with that of an actual missile. The measured time constant 

of the missile was approximately 0.21 seconds over the various 

flight conditions. The measured response of the model was 

0.19 seconds over the same flight conditions. 

The transient response time to a step torque, equivalent 

to a step input (<5 ) of 3° was approximately 0.5 seconds for 

c 

the actual missile and 0.62 seconds for the model (shown in 
Figure 2-5 ) . 

Tests were conducted on the missile model under two sets 

of flight conditions covering the limits of normal operation. 

Flight condition No. 15 represented normal altitude and 

velocity, whereas flight condition No. 24 represented low 

altitude and normal velocity. The self-adaptive nature of 

the control loop can be seen in Figures 2-4 through 2-10. 

The amplitude of K^F^ changes according to the aerodynamic 

conditions, but the response (<p/<p or (p/6 ) remains constant. 

c c 
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FIGURE 2-2 



($>' vs TIME with No Inputs 




Xscale: 1.0 sec/inch 
Yscale: 4°/inch 



Note: Initial excursion is due to the initial conditions 
used in the simulation. 
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FIGURE 2-3 



<j> r vs TIME with -45° step input @ t=3.0 sec 
( flight condition #15) 




2.3: j: a. a: 3.2: 12.:: 



Xscale: 2.0 sec/inch 
Yscale: 10°/inch 
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Curve #2 Step Input 
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FIGURE 2-4 



K F vs TIME with -45° step input @ t=3.0 sec. 
(^f fight condition #15 ) 




Xscale: 2.0 sec/inch 
Yscale: 20°/sec/inch 
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FIGURE 2-5 



$1 vs TIME with +3° step torque (<5 c ) at t-3.0 sec. 
( flight condition #15 ) 



rt c* 




Xscale: 2.0 sec/inch 
Yscale: 4°/inch 
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FIGURE 2-6 



K F vs TIME with +3° step torque (6 ) at t=3.0 sec. 
( 6 f fight condition #15 ) 




Xscale: 2.0 sec/inch 
Yscale: 40°/sec/inch 
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FIGURE 2-7 



({('vs TIME with a -45° step input at t-3.0 sec. 
( flight condition #24 ) 




Xscale : 2.0 sec/inch 
Yscale: 10°/inch 



Curve #1 <{>' 

Curve #2 Step Input 
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FIGURE 2-8 



K F vs TIME with a -45° step input at t=3.0 sec. 
(“flight condition #24 ) 




Xscale: 2.0 sec/inch 
Yscale: 40°/sec/inch 
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FIGURE 2-9 



<j)'vs TIME with a +3° step torque (6 ) at t=3.0 sec. 
( flight condition #24 ) 




Xscale: 2.0 sec/inch 
Yscale: 10°/inch 
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FIGURE 2-10 



K F vs TIME with a +3° step torque (6 ) at t=3.0 sec. 
(^flight condition #24 ) 




Xscale: 2.0 sec/inch 
Yscale: 100 °/sec/inch 
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D. CONCLUSIONS 



The simulation flights shown in Figure 2-2 through 
Figure 2-10 were compared to the results of test flights 
performed by the civilian contractor. On the basis of these 
comparisons, the conclusion was reached that the model was an 
accurate representation of the actual missile control system 
and that the model could therefore be used in the study of 
the roll wander phenomenon. 



III. ROLL WANDER 



A. INTRODUCTION 

As stated earlier, the two major areas of research were 
the construction and testing of an accurate computer simula- 
tion for the missile and the study of the phenomenon of "roll 
wander. " 

This undesired roll loop excitation or roll wander re- 
sults from noise originating in the tail control surface 
actuators. This noise is sensed by the roll rate gyro and 
is thus introduced by the roll rate gyro into the rate feed- 
back loop as a roll rate error signal. The production of 
roll wander then occurs in one of two fashions. 

B. LOW AMPLITUDE ROLL WANDER 

The noise signal introduced is enhanced by the effects 
of the roll loop nonlinearities. The effect of these non- 
linearities is to generate harmonics of the dither frequency. 
These harmonics in turn add with the fundamental dither fre- 
quency to produce a distorted wave form of a frequency equal 
to the fundamental frequency but with zero crossings within 
one cycle unevenly spaced (Figure 3-1) • Since the bistable 
element is triggered only by these zero crossings, the average 
frequency of the bistable remains the same; but the time in 
its two states becomes unbalanced (Figure 3-2). The neutral 
position of the tail surfaces is the integration of the 
bistable element output. Therefore, the tail control surfaces 
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FIGURE 3-1 



K F vs TIME Showing Harmonic Addition 
(“Harmonic commences at t=1.0 sec. ) 




Xscale : 0.02 sec/inch 
Yscale: 20°/sec/inch 



#1 Resultant 'distorted waveform 
tr 2 Fundamental dither frequency 
#3 Harmonic frequency 
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FIGURE 3-2 



Bistable Output 
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begin to move off of their zero positions and the missile 
rolls. Since no roll command was introduced into the loop, 
a resultant dc error signal is generated in the roll position 
loop (outer loop). This dc signal causes the input to the 
bistable element to be shifted accordingly above or below 
the zero axis. Since the distortion of the waveform is 
caused by harmonic addition, and the fundamental frequency 
of the bistable remains unchanged, a shifted position exists 
for the waveform where the zero crossings once more become 
equal. The result is that the missile initially rolls off 
zero upon introduction of the harmonic, then rolls back 
(Figure 3-3) • If the phase of the harmonic wave is shifted 
rapidly enough to repeatedly offset the resultant waveform, 
the missile will begin to wander in roll at a frequency of 
1.2 to 2.5 hz trying to compensate (Figure 3 _i *)- In the 
actual missile, uncorrelated white noise from the tail con- 
trol actuators produces this phenomenon. The broad band 
noise contains components at the harmonic frequencies of the 
dither signal. Since this noise is random uncorrelated 
noise, there exists no coherence in the phase of the various 
spectral components. Therefore, the action of the noise 
is the same as the addition of the fundamental dither fre- 
quency with a series of harmonic frequencies shifting 
rapidly in phase. Noise with very low amplitude, only 10 % 
of the amplitude of the signal (SN = +20db) is sufficient 
to cause the onset of roll wander (Figure 3-5) • 



FIGURE .3-3 



<j>' vs TIME Showing the Effect of a Fixed Phase Harmonic 
( Harmonic Commences at t=1.0 sec. ) 







Xscale: 1.0 sec/inch 
Yscale: 4°/inch 



Note: Initial excursion is due to the initial conditions 
used in the simulation. 
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FIGURE 3-4 



<p' vs TIME Showing the Effect of a Varying Phase Harmonic 
( Harmonic Commences at t=1.0 sec. ) 




Xscale: 1.0 sec/inch 
Yscale: 4°/inch 
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FIGURE 3-5 



4 / vs TIME Showing the Effect of Adding White Noise 
( constant level ) Noise Commences at t=0.7 sec. 




Xscale: 2.0sec/inch 
Yscale: l°/inch 

Note: Initial excursion is due to the initial conditions 
used in the simulation. 



C. LARGE AMPLITUDE ROLL WANDER 



As the noise amplitude increases, the amplitude of the 
wander increases, since the harmonic components of the noise 
signal also increase. The frequency of the wander, however, 
remains constant (Figure 3-6). The action of the noise is 
to reduce the low frequency roll loop gain. The decrease 
in the low frequency roll loop gain will decrease the roll 
system damping in the closed loop response, resulting in 
larger wander excursions. If the noise gets large enough, 
the bistable switch gain may be decreased to where the roll 
system will enter into a sustained low frequency oscillation 
of extremely large magnitude. The dithering action of the 
loop is. stopped and the entire system operates in a saturated 
limit cycle mode of operation (Figure 3-7) • 

D. SIMULATING THE ROLL WANDER 
1 . Noise Source 

One of the more difficult problems was that of accu- 
rately simulating the roll wander phenomenon. Analog simula- 
tions performed by the civilian contractor had been performed 
using *100 hz band limited gaussian noise. To produce this 
digitally, a FORTRAN subroutine was written utilizing a 
random number generator and a digitally simulated *100 hz low 
pass filter to generate the appropriate band limited noise. 
The digital simulation of the transfer function of a low 
pass filter requires the integration of all data points to 
produce the filtered output. Therefore, to preclude the 
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FIGURE 3-6 



<f>' vs TIME Showing the Effect of Adding White Noise 
( linearly increasing level ) Noise Commences at t=0.7 sec 




Xscale: 2.0 sec/inch 
Yscale: 5°/inch 
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FIGURE 3-7 



4>* vs TIME Showing Loop Instablity due to Saturating 
the Loop with High level White Noise 



„■ 1 * t 

C» C? 




Xscale: 2.0sec/inch 
Yscale: 2.0x10 “/inch 

Note: Roll loop enters saturation at approximately 
t=8.0 sec. 
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necessity of those time consuming integrations during the 
actual missile simulation, twenty thousand points of 400 hz 
bandlimited white noise were computed, then written on disk. 

A subroutine in the DSL program called NOISE was called each 
time increment. With each call to subroutine NOISE a read 
was performed from disk and the required white noise w as 
supplied as output. Since a DELT (integration step size) of 
0.0001 was used in the simulation, two seconds of white 
noise were available. At the end of that time, the noise 
record was rewound and reused. The effect of this repetition 
of the noise record every two seconds of simulation time on 
the randomness of the white noise was considered negligible. 

A sample of the noise used in producing the wander shown in 
Figure 3-6 is shown in Figure 3-8. 

2 . Harmonic Generator 

Since low amplitude wander vice large amplitude roll 
wander was of prime importance, it was decided, that to 
effectively study the phenomenon with the hope of preventing 
it, a pure harmonic waveform rather than broad band noise 
should be used as an input. 

An "in-phase" second harmonic could be generated by 
using a square law device . 

HARM = BI 2 

then if: BI = BI COS 9 

max 
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FIGURE 3-8 



White Noise (linearly increasing level ) 




). i :: s. j; t. r, jj. : 



Xscale: 2.0 sec/inch 
Yscale: 80°/sec/inch 
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HARM = BI 2 COS 2 0 



BI 2 (?+4Co S-6 ) 

max 2. d 



max 



Bf = BI + HA RM 



Bf = BI 



max 




( a + a Cos?e ) 



(15) 



The harmonic term however carries with it a dc term which 
would cause a roll offset, and more importantly, might hide 
the dc error term generated by the control loop. To eliminate 
this dc level, the output HARM1 was passed through a three 
stage low pass filter with a -3db frequency of 10 hz and 
a -18 db/oct roll off. 



offset exists, a dc error term is produced to null out the 
loop. When the harmonic is added to the loop, an immediate 
dc term is produced moving the fundamental waveform off the 
null axis. Since 'the harmonic is obtained by squaring the 
signal in the outer loop, if this signal has such a dc error 
term then: 



As discussed earlier in the chapter, when a roll 



HARM = BI 



2 



(BI + BI X COS 0) 



2 



BI + 2BI -BI, COS 0 
o o 1 



(16) 



+ BI 2 (^ + ^ Cos a©) 
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Equation (16) represents a distorted waveform. As 
BI (the dc component) increases, the effect of the harmonic 
component decreases. To avoid this problem and again to 
enable viewing only the wander phenomenon, the waveform used 
to produce the harmonic was passed through another three 
stage low pass filter similar to the one described above, 
(Figure 3-9). Thus at anytime, only a pure sinsoid at a 
harmonic frequency was introduced. 

Addition of this phase locked harmonic is shown in 
Figure 3-2. The excursion caused by the addition of a har- 
monic signal producing a signal to noise (harmonic) ratio 
of approximately +12 db was -11.18°. To produce sustained 
roll wander, the phase of the harmonic signal had to be 
shifted rapidly in relation to the fundamental. To accom- 
plish this, a DSL DELAY function block was used. The amount 
of the delay was varied by setting it equal to: 

Shift = r Sih S TT 

TO 

where: ^ 

4 S = 10 hz 

4o = fundamental dither frequency 
then: HARM’ ( s ) = HARM (S) e ~^ shift ^ s 

The resulting waveform is shown in Figure 3-3- The phase 
shifting frequency used in this simulation was fast enough 
to cause an increasing amplitude wander to be produced. 



SIMULATION LOOP with HA.RHONIC GENERATOR 
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E. CONCLUSIONS 



Producing the roll wander with random noise signal 
ensured that the model was behaving properly. The model 
reacted to the introduction of noise in the same fashion 
as the real missile. 

The use of the harmonic generator to produce a sustained 
smooth roll wander proved that the phenomenon was in fact 
caused by harmonic addition. In addition, the use of the 
harmonic generator gave valuable insight into the roll 
loop’s response to external stimulus. The fact that the 
output of the bistable element drives the control surfaces 
at such a high rate, producing such large missile roll rates, 
explains why it is so difficult to view an imbalance in 
the bistable output. 
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IV. A SOLUTION 



A. INTRODUCTION 

Once the cause of the roll wander phenomenon had been 
isolated and verified through testing of the missile model 
with the noise source and harmonic generator, the problem 
existed to find a way of preventing this wander. As dis- 
cussed earlier, two causes of wander due to a noise input 
exist. Measures to prevent large amplitude wander were not 
considered. Large amplitude wander would be encountered 
only under extremely abnormal conditions in which the noise 
signal totally saturates the control loop. On the other hand, 
low amplitude wander caused by harmonics present in the noise 
adding with the fundamental dither frequency, was felt to 
be a significant problem. In this chapter, the term roll 
wander will be understood to deal with low amplitude wander. 

The most direct method of preventing roll wander would 
be the attenuation of the dither harmonics in the roll con- 
trol loop. This was a deceptively difficult problem. There 
already existed in the control loop several high frequency 
filters to filter out those high frequency signals generated 
by the missile's elastic body bending modes. Notch filters 
vice a simple low pass filter had been utilized to prevent 
the shifting of the fundamental dither frequency due to the 
phase shift imposed by the addition of the filters. Inclu- 
sion of another set of filters notched at the main harmonic 
components of the dither fundamental could be used to combat 
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the frequency addition and resulting roll wander. However, 
since the dithering action of the roll loop is dependent 
upon the phase crossover of the loop, which is determined by 
the loop elements, the frequency of the dither is dependent 
on the roll loop parameters and varies from missile to 
missile. More importantly, the dither frequency varies 
slightly from one set of flight conditions to another. Thus, 
a notch filter set to filter out harmonics would of necessity 
have to have a notch wide enough for all frequency variations. 
This, in turn, would cause an unacceptable decrease in the 
dither frequency. A self-adaptive notch filter was thus 
proposed as a possible solution. 

B. NOTCH FILTER DESIGN 

A filter with conjugate imaginary zeroes located at the 
frequency of interest, the notch frequency, would provide 
total attenuation of the signal at that frequency. At fre- 
quencies higher than the notch frequency, the filter would 
have a gain greater than unity. However, the filter must 
have unity gain at all frequencies other than the notch 
frequency. To counter the gain of the filter produced by 
the conjugate imaginary zeroes, a double pole is required. 

For stability reasons, purely imaginary poles are not possi- 
ble. The poles must have an attendant 5 term. A problem 
arises, however, when complex conjugate poles with non zero 
real parts are introduced to produce a realizable filter. 

There is a sign reversal at the notch frequency, which 
implies a l80° phase shift about the notch. Since the zeroes 
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are pure imaginaries, having no real parts, there is no 
phase disturbance at other frequencies due to this shift. 
This is not the case with the complex poles. The value of 
C in the denominator of the filter expression determines the 
amount of phase disturbance the filter produces . Frequency 
components at lower frequencies, specifically at the funda- 
mental frequency, receive a phase shift depending on the 
magnitude of £ . The shift in phase produced by the filter 
changes the phase crossover point of the entire loop. The 
phase crossover sets the limit cycle oscillation (dither) 
frequency. The value of £ also determines the effective 
width of the filter notch. An engineering trade-off is 
required, therefore, between the width of the notch and the 
change produced in the dither frequency. 

Finally, a pole must be added at some higher frequency 
(3 = 3000) to allow the response of the filter to roll off 
at high frequencies. 

The judgment was made in the design of this filter that 
the zeroes of the filter were to be pure imaginaries. This 
does not imply that complex conjugate zeroes with non zero 
real parts would not also provide attenuation at the notch 
frequency. However, if the zeroes also had non zero real 
components, the amount of attenuation about the notch fre- 
quency would be decreased. With the width of the notch 
already a critical factor; and since imaginary zeroes are 
possible (in practice; zeroes with negligible real parts), 
the filter was designed as stated. Equation (16) is the 
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form of the transfer function of the filter proposed. 



S* + (Jo 2 (17) 

{ S 2 4 2 2 S + ) ( 5 -v (S ) 

Figures 4-1 through 4-6 show the Bode magnitude and phase 
plots for the filter's frequency response from 1 hz to 1.59 
khz for values of c; of 0.01,, 0.1 and 0.2 respectively. The 
magnitude plots show the notch formed at the center frequen- 
cy. (A frequency of 56 hz, twice the fundamental dither 
frequency was used.) From the phase plots can be seen the 
amount of phase shift introduced by the filter. With a C 
of 0.2, a phase shift of -25° was experienced at 28 hz. 

For £ equal to 0.1, the phase shift was -17.7°; and with C 
equal to 0.01, the phase shift was -10.8°. 

The filter was tested outside of the missile model. It 
was found that adequate attenuation of a signal having a 
bandwidth of 0.2 hz centered at the filter notch frequency 
was not obtained with values of £ less than 0.07. When the 
filter was tested in the missile model with C set to 0.01, 
the frequency of the dither was reduced from 28.75 hz to 
26.3 hz; with £ set to 0.1, the frequency of the dither 
dropped to 24.5 hz, and with a C of 0.2, the frequency 
dropped to approximately 22 hz. 

It was decided to use a value for t; of 0.1 for minimum 
disturbance in the loop with maximum attenuation of the 
harmonic frequencies. Two filters were inserted, in cascade. 
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FIGURE 4-1 



Bode Magnitude Plot of the Notch Filter with £--0.01 




Xscale: log ( lOhz ) /inch 
Yscale: 8.0x10 3 • units/inch 
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FIGURE 4-2 



Bode Phase Plot of the Notch Filter with c=0.01 




Xscalo : log ( IChz ) /inch 
Yscale: 40°/inch 
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FIGURE 4-3 



Bode Magnitude Plot of the Notch Filter with £=0.1 
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FIGURE 4-4 



Bode Phase Plot of the Notch Filter with £=0.1 



<> n 

t l 




Xscale: log ( lOhz ) /inch 
Yscale: 40°/inch 
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FIGURE 4-5 



Bode Magnitude Plot of the Notch Filter with r,- 0.2 




Xscale: log ( 10h? ) /inch 
Yscale: 8.0x10 0 units/inch 
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FIGURE 4-6 



Bode Phase Plot of the Notch Filter with £=0.2 




Xscale: log ( lOhz ) /inch 
Yscale : 40°/inch 



66 



with notch frequencies set to filter out the second and 
third harmonic frequencies of the fundamental. 

C. SELF ADAPTIVE FEATURE 

The notch frequency of the filter is specified by the 
value of u> o in Equation ( 16 ). The compromise concerning the 
value of c, which determines the width of the notch was 
based upon the assumption that would be within ±0.1 hz 
of the harmonic frequency to be filtered. However, as 
previously discussed, several factors, including the inclu- 
sion of the filter itself into the loop, will cause the 
fundamental frequency to vary. When this happens, the 
harmonic of the fundamental also varies, by a larger amount, 
proportionate to the order of the harmonic in question. If 
the harmonic component moves in frequency outside the .2 hz 
wide notch, the filter no longer has any effect. The filter 
must therefore be self-adaptive; it must have a variable 
notch frequency. 

The filter obtains its self-adaptive feature by first 
measuring the dither frequency, then setting the filter 
coefficients to produce the correct notch. The filter pro- 
duces an initial transient response. Constantly changing 
the filter coefficients, which is effectively the same as 
constantly inserting a different filter, disturbs the control 
loop. This disturbance is seen as a random wander in the 
position loop. To minimize this effect, the filter coeffi- 
cients must be held constant or be required to change slowly 
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in a smooth fashion. This will only be the case if the 
dither remains fairly stable. 

A jitter in the dither frequency is produced in the 
simulation because continuous integrations are being simu- 
lated by discrete numerical operations. Also, since the 
output of the bistable element is a pure square wave , 
instantaneous frequency was most easily measured at this 
point. Any external influence which disturbs the relative 
sizes of the states of the bistable element would also 
produce a fluctuation in the instantaneously measured 
frequency. The measured value of instantaneous frequency 
did not therefore remain constant. When the filter coeffi- 
cients were changed each integration interval, to correspond 
to this instantaneous frequency, a wander with a frequency 
of approximately 2.5 hz was observed. The filter was itself 
causing a roll wander phenomenon (Figure 4-7). 

This side effect was significantly reduced by using an 
average value for the frequency. The frequency was measured 
each integration interval, then averaged over six cycles. 

As each new cycle was recorded, it was then averaged with 
the previous five cycles; thereby "fairing" out any frequency 
jitter. The result is seen in Figure 4-8. 

It is obvious that in actual practice the construction of 
such a self-adaptive filter is not quite so simple as the 
simulation of one. However, techniques do exist for the 
production of such filters, either analog or digital. This 
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FIGURE 4-7 



<j>'vs TIME with Self Adaptive Filter in Control Loop (no input) 
Adaptation using Instaneous Frequency 




Xscale : 2.0 sec/inch 
Yscale: 2.0°/inch 

Note: Initial excursion is due to the initial conditions 
used in the simulation. 

Note: Dither frequency is measured each integration interval. 
Filter coefficients are updated at the same rate. 
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FIGURE 4-8 



<p' vs TIME with Self Adaptive Filter in Control (no input) 
Adaptation using Averaged Frequency 




13 . , 







Xscale: 1.0 sec/inch 
Yscale: 4.0 “/inch 



Note: Dither frequency is measured and averaged over 6 
cycles. Filter coefficients are updated using averaged 
value of frequency. 
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thesis concerned itself only with the filter’s effect on 
the control loop, and not on filter construction. 

D. TESTING THE FILTER 

Since it was the effect of the addition of harmonic 
frequencies with the fundamental dither frequency which 
produced the distorted waveform responsible for roll wander, 
both the harmonic generator and the filter were inserted 
into the loop just prior to the bistable element (Figure 
4-9). In the actual missile, it would be expected that such 
a filter would be placed at the output of the roll rate gyro; 
since it is the roll rate gyro which inputs the actual 
harmonic frequencies into the loop. 

With the filter installed, the control loop model was 
tested in the same fashion as described in Sections II and 
III. Figure 4-8 and Figures 4-10 through 4-l4 show the 
response of the model to no input, step input and a tail 
induced torque, respectively. These runs were all made with 
normal, identical- flight conditions. The missile’s response 
during high altitudes and speeds remained the same as before 
inclusion of the filter. 

E. WANDER SUPPRESSION 

The runs just described showed that with the self- 
adaptive filter in the loop, the model still responded 
properly to input commands. The model was then tested again 
with the harmonic generator. Figure 4-15 and Figure 4-l6 
clearly show the filter’s effect. The roll wander is totally 
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FIGURE 4-10 



K F. vs TIME with Filter in Loop (no input) 
6 5 




Xscale: 2.0 sec/inch 
Yscale: 20°/sec/inch 



73 






FIGURE 4-11 



<f> ‘ vs TIME with Filter in Loop and -45° step input at t=3.0 sec 




Xscale: 2.0 sec/inch 
Yscale: 10°/inch 

Curve #1 <j>' 

Curve #2 Step Input 
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FIGURE 4-12 



K^F~ vs TIME with Filter in Loop and a -45° step input 
at t=3.0 sec. 




i.y» t. j: 3.3: jj. :: 



Xscale: 2.0 sec/inch 
Yscale: 20°/sec/inch 
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FIGURE 4-13 



<*/ vs TIME with Filter in the Loop and a +3° step torque (6 ) 
at t=3.0 sec c 



U O 




Xscale: 2.0 sec/inch 
Yscale: 4°/inch 
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FIGURE 4-14 
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Xscale: 2.0 sec/inch 
Yscale: 40°/sec/inch 
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FIGURE 4-15 



<$' vs TIME Showing the Effect of a Varying Phase Harmonic 
( Harmonic Commences at t=1.0 sec ) 



C 1 L? 

p i >4 




Xscale : 1.0 sec/inch 
Yscale : 4° /inch 

Note: Initial excursion is due to the initial conditions 
used in the simulation. 
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FIGURE 4-16 



</ vs TIME with Filter in Loop and with a Varying Phase 
Harmonic ( Harmonic Commences at t=1.0 sec. ) 




Xscale : 1.0 sec/inch 
Yscale: 4°/inch 

Note: Initial excursion is due to the initial conditions 
used in the simulation. 

Note: The excursion at t=1.0 sec. is due to the harmonic 
causing a momentary shifting of the dither frequency outsid.e 
the notch bandwidth. 
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suppressed by the notch filter. The transient disturbance 
seen in Figure 4-l6 is caused by the initial "switching in" 
of the harmonic signal at t equals one second. The initial 
effect of the harmonic produces a change in the dither 
frequency, with an attendant change in the output signal of 
the harmonic generator. Since the notch filter is using an 
averaged value of frequency in computing the filter coeffi- 
cients, the shift in harmonic frequency is momentarily 
outside the range of the notch. Roll wander starts. As 
the averaged frequency approaches the actual loop frequency 
the notch shifts and the wander is suppressed by the filter. 

Figures 4-17 and 4-l8 show the model's response to the 
noise input. A noise signal will disturb the roll loop and 
cause the missile to roll. This, however, is not roll 
wander. Comparing these two figures, it can be seen that 
the amplitude of the wander has been decreased by about 
50 per cent. The rougher appearance of the curve in Figure 
4-l8 is due to the fact that the random disturbance of the 
noise is becoming more dominant as the harmonic addition is 
being suppressed. 

F. CONCLUSIONS 

The primary objectives of this thesis were: 

1. Construct an accurate digital model of the missile 
roll control loop. 

2. Produce in the model the phenomenon of roll wander 
with an explanation of the sources and mechanisms of 
the wander. 

3. To produce a means of preventing or at least signi- 
ficantly reducing the occurrence of this phenomenon. 
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FIGURE 4-17 



^ vs TIME Showing the Effect of Adding White Noise (constant 
level). Noise Commences at t=0 . 7 sec. 




Xscale: 2.0 sec/inch 
Yscale: l°/inch 

Note: Initial excursion is due to the initial conditions 
used in the simulation. 
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FIGURE 4-18 



<t>' vs TIME with Filter in Loop and White Noise Added 
Noise Commences at t=0.7 sec. 




Xscale: 2.0 sec/inch 
Yscale: l°/inch 
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C. 



The first two objectives have been achieved in this 
thesis. As shown in Sections I, II, and III, the model 
behaves in the same fashion as the real missile. The digital 
simulation of the missile now exists at the Naval Postgraduate 
School, and is available for further study. 

The solution presented in this chapter as a cure to the 
roll wander problem works and will prevent this phenomenon. 
However, a self-adaptive filter has several weaknesses which 
limit its effectiveness. The inclusion of the filter into 
the control roll has reduced the fundamental dither frequency 
significantly. Such a reduction may cause interference 
problems with the missile guidance system. It was found that 
changing the filter coefficients, to change the notch fre- 
quency, disturbed the control loop, causing a fluctuation in 
the fundamental frequency and random wander in the position 
loop. To minimize this effect, an averaged frequency was 
used. Because of this, the filter did not respond to fre- 
quency changes as quickly; and the roll wander caused by the 
noise input was not totally suppressed. 

It is believed that much room exists for future thesis 
work in the area of preventing the roll wander. The largest 
portion of the work involved in this thesis was the construc- 
tion of an accurate model and the production of the roll 
wander. With the model now in existence, future thesis 
students can use it and concentrate on the development of 
a positive efficient system to prevent roll wander. 
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V. RECOMMENDATIONS FOR FUTURE STUDY 



As discussed in the previous chapter, continued study 
is needed in the implementation of an adaptive filter for 
the suppression of roll wander. The filter used in this 
thesis produced a decrease in dither frequency of almost 
4 hz. This alone may make the filter acceptable. 

Another area of study, apart from the study of the missile 
control system, also became evident during the course of this 
thesis. That is the study and the use of digital simulation 
languages. As is often the case with any research, the 
preparation consumed as much time as the investigation itself. 
I am referring to the production of the DSL plot package. 

DSL lay dormant at the Naval Postgraduate School because of 
its inability to output to the CALCOMP plotter. Once an 
operational plot package was produced, the value and the 
power of DSL as a simulation tool rapidly became apparent. 

IBM’s Continuous System Modeling Program (CSMP), another 
simulation language at this school, also at present requires 
considerable programming to access the CALCOMP plotter. It 
would require only slight modification to DSLPLOT and minor 
procedural changes to the system control cards for CSMP to 
make CSMP output to the plotter via DSLPLOT. 

It must be remembered that as with any simulation routine 
DSL is time consuming; and on a time-shared system such as 
the IBM 360/67 System used at the Postgraduate School, turn 
around time is, at best, marginal. The IBM 360/67 System 



has interactive graphics terminals. If DSL were modified 
via appending FORTRAN programs to output and input via a 
graphics terminal, a pseudo real time system could be 
produced. The advantages of such a system would be great. 

The student could modify his system, and automatically see 
the result. If the program misbehaved, the student could 
terminate the run immediately, rather than waiting the full 
specified time of the run and the time required to exit the 
system prior to being able to study the program and take 
corrective action. Such a system could easily allow the 
user to accomplish in a two or three hour session what might 
require days or weeks to accomplish off line. The modifica- 
tion of DSL to produce such a system would be extensive but 
rewarding. 

The other modeling languages on file at the Postgraduate 
School (CSMP, ECAP, LISA, etc.) all would benefit immensely 
by similar modifications. The saving in systems time by 
removing such jobs from the normal job queues and rescheduling 
them during hours of minimum load would also be significant. 

The emphasis in control engineering has shifted dramati- 
cally to the utilization of the digital computer, both for 
problem solving and system modeling. However, the structure 
and implementation of the simulation programs available is 
not explicitly taught. A course built around system modeling 
using such programs would be highly beneficial. 
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APPENDIX A 



A . INTRODUCTION 

Digital Simulation Language (DSL) is an IBM System/360 
FORTRAN IV program for the digital simulation of continuous 
system dynamics. Its nonprocedural problem-oriented input 
language accepts problems expressed either at the analog 
block diagram level or as systems of ordinary differential 
equations. Since the input language is nonprocedural, proper 
statement sequence is established by an internal sort based 
on correct information flow. DSL includes FORTRAN as a 
subset, thereby extending its power to handle nonlinear and 
time-variant problems of considerable complexity. A trans- 
lator converts the input statements into a FORTRAN subprogram, 
which is compiled and executed to accomplish the simulation. 

A significant feature of DSL is centralized integration, 
which means that all integrator outputs are computed simul- 
taneously at the end of each calculation interval. The user 
may specify one of several integration methods, including 
fifth order Milne Predictor and Corrector and second and 
fourth order Runga Kutta. 

DSL as originally presented at the Naval Postgraduate 
School, outputs to the high speed printer only, either as 
printed output or printer plots . Due to software incompat- 
ibilities, DSL as initially installed would not output to 
the CALCOMP XY plotter. This thesis, as with many simulation 
studies, required the CALCOMP as an output. For this reason. 
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and since DSL had so many advantages, extensive work was 
done prior to constructing the missile simulation on a FORTRAN 
plotting program to output the DSL program via the CALCOMP 
plotter. This program, DSLPLOT and its associated FORTRAN 
program PLOTHALL, will be discussed later. 

B. PROGRAM FLOW 

At the start of the program, all program variables are 
set to the initial conditions given. TIME, the independent 
variable, is set to zero, and a run counter KSIM is incre- 
mented if it is the beginning of a new parameter study. 
Initially KSIM is set equal to 1. At this point, the region 
of the program called the INITIAL Region is entered. In the 
INITIAL Region, arithmetic operations and logical statements 
to be executed only at the beginning of a particular run are 
computed. An integer flag KEEP, is set equal to zero, and 
the DERIVATIVE Region is called. 

The DERIVATIVE Region represents that portion of the 
simulation which involves integration and the calculation of 
the derivatives. The basic interval in the independent 
variable for each pass through this region is the calculation 
interval DELT. It is for this reason that the DERIVATIVE 
Region should include all time (independent variable) 
dependent calculations which must be performed at the same 
frequency required by the selected integration algorithm to 
assure accurate numerical integration. 
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The nonprocedural nature of the input language is the 
result of the DERIVATIVE Region. The DERIVATIVE Region is 
the only region in which an internal sort routine is utilized 
to establish the correct statement sequence. Only structure 
statements are sorted and sequenced. A statement is properly 
sequenced if all its inputs are available either as input 
parameters or initial conditions or as previously computed 
values in the current calculation interval. Outputs of INTGRL 
and MEMORY blocks are known because their initial values are 
given. The internal sequencing procedure begins first with 
integrator outputs, then with outputs of MEMORY blocks and 
finally with those remaining output variables in the 
DERIVATIVE Region that have not been included in the sort 
sequences of integrators and MEMORY blocks. It is important 
to realize that the sorting algorithm is unable to sequence 
the statements properly unless it can assume a known output 
variable as a starting point of a sort sequence. As previously 
stated, this output is provided by either an INTGRL or MEMORY 
block. A loop which does not contain at least one such block 
is called an IMPLICIT or ALGEBRAIC loop and cannot be 
sequenced. Upon encountering such a loop, the sorting algo- 
rithm is aborted, and all subsequent statements are executed 
in order. 

The use of the sort algorithm places a severe restriction 
on the type of statements which can be used in the DERIVATIVE 
Region. No FORTRAN logic or branch statements can be used 
because they cannot be sorted properly. To allow the 
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versatility which accompanies these types of statements, 
once every integration interval, the program flow exits the 
DERIVATIVE Region and enters the DYNAMIC Region. The DYNAMIC 
Region provides an area in which FORTRAN logic can be used. 
Testing of variables, adjustment of coefficients and nonlinear 
operation all can be accomplished. It should be remembered, 
however, that since for a given value of the variable TIME, 
the DERIVATIVE Region is executed before the statements in 
the DYNAMIC Region a time lag equal to DELT exists between 
the actions taken in the DYNAMIC Region and the resulting 
effect in the DERIVATIVE Region. 

A second system variable, DELS, dictates the frequency 
at which the computation interval is interrupted and program 
flow is diverted to the SAMPLE Region. The makeup of the 
SAMPLE Region is the same as that of the DYNAMIC Region. 

Its main function is to test variables to determine points 
for data collection and to order the program flow to allow 
simultaneous collection of data of several variables. 

DSLPLOT, which was installed at the Naval Postgraduate School, 
is accessed via the SAMPLE Region. 

Finally, at the termination of the simulation either by a 
logic statement or when the variable TIME is set equal to the 
system variable FINTIM (a variable set by the user) , the 
program makes one last pass through the SAMPLE Region and 
enters the TERMINAL Region. In the TERMINAL Region, logic 
can be used to modify program variables and constants and a 
new run initiated. If DSLPLOT has been utilised, the plotting 
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data set is closed out. If no action is taken to initiate 
a new run, the DSL job will be terminated. If a new run 
has been called for program flow returns to the INITIAL 
Region . 

When the DSL job is completed, the DSL program is 

dynamically linked via the IBM 360 linkage editor to the 

program DSLPLOT. 

/ 

C . DSLPLOT 

The program DSLPLOT was produced as a forerunner to this 
thesis to allow the DSL output variables to be displayed via 
the CALCOMP XY Plotter. DSLPLOT accomplishes this by writing 
the variables to be plotted along with counters indicating 
the plot number, curve number, and individual point number on 
an external data set via the subroutine DRWG each sample 
interval (DELS). Each group of data points indicating a 
single DSL run is delineated by a set of zeroes produced by 
the subroutine ENDRW called in the TERMINAL Region. 

PLOTHALL, the- appending FORTRAN program used by DSLPLOT, 
accesses this data set upon a successful completion of the 
DSL job. A three-way sort is performed on the data set, and 
the points to be plotted are written in order on a second 
data set. The points are then retrieved one curve at a time 
from the second data set and passed to the plotting subroutine. 
This subroutine utilizes the IBM Plotting Package and produces 
the required CALCOMP plots. 
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Since it was intended that DSLPLOT was to be a permanent 
addition to the DSL program at the Naval Postgraduate School, 
it was written with as much flexibility as possible. Listed 
in Table A-l are the various modes under which DSLPLOT can 
operate. The different modes allow the user to choose the 
size of the output plot and to either specify the scaling 
used on either axis or to allow the program to auto scale one 
or both axes . 

D. USE OF THE PLOT PACKAGE 

As previously discussed, points are collected for plotting 
via the subroutine DRWG. A call to this subroutine is 
initiated by the FORTRAN statement: 

CALL DRWG (1,1, NUM,X,Y) 

Where X and Y are the variables to be plotted. The X 
variable is plotted as the abscissa coordinate and Y the 
ordinate value . MUM is the integer counter denoting the data 
point number. The two integer constants appearing as the 
first two arguments represent the plot number and curve 
number for that plot number. These arguments may be provided 
as integer variables. This allows overlaying a set of curves 
from different runs on the same plot. It must be remembered 
that there is no convention in DSL, as there is in standard 
FORTRAN, for providing integer variables. Those variables 
such as NUM, or the variable NPLOT used in ENDRW must be 
declared integers by a DSL INTGER statement. 
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TABLE A-l 



Mode No . User Supply (in order listed) 



0 


Auto 


(no 


inputs 


required) 




1 


Xmin 


Dx 


Ymin 


Dy 






2 


Xmin 


Dx 










3 


Ymin 


Dy 




















Lenx 


Leny 


5 


Xmin 


Dx 


Ymin 


Dy 


Lenx 


Leny 


6 


Xmin 


Dx 






Lenx 


Leny 


7 


Ymin 


Dy 






Lenx 


Leny 


8 










Lenx 




9 


Xmin 


Dx 


Ymin 


Dy 


Lenx 




10 


Xmin 


Dx 






Lenx 




11 


Ymin 


Dy 






Lenx 




12 










Leny 




13 


Xmin 


Dx 






Leny 




Ik 


- Xmin 


Dx 


Ymin 


Dy 


Leny 




15 


Ymin 


Dy 






Leny 





Note: If mode 0 is specified a 9" by 9" graph with all 

values auto scales will be produced. 
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Upon completion of data collection for a given set of 
plots, the data set is closed out by the FORTRAN statement: 

CALL EN DRW (N PLOT) 

Where NPLOT is an integer variable which is initially set to 
the number of intended data sets. Each call to ENDRW 
decrements the variable by one. When NPLOT is less than or 
equal to zero, the DSL job is terminated. After a call to 
ENDRW, if NPLOT is still greater than zero, the program 
simply exits the subroutine, and the DSL program continues. 

Each data set must be closed out by a call to ENDRW. 

If the plotting package has been used in the job, ENDRW must 
also be used as the means of terminating the DSL job. 
PLOTHALL, link edited to the DSL job, terminates upon encoun- 
tering a value of NPLOT equal to one, therefore if the main 
DSL job is terminated by means other than ENDRW, an abnormal 
termination of the plot package will result. 

Since it is standard practice to initiate a new set of 
plots each DSL run, the plot, curve and data point indicators 
refer to a given run. After a call to ENDRW, if a new run is 
initiated either as a parameter study or by a call to RERUN, 
CONTIN, etc., encountering the same DRWG statement again 
will generate a new plot. If curves are to be overlaid on 
one plot; at the end of a given run the data set must be held 
open by not calling ENDRW and the curve indicator, specified 
by a variable name, is updated. When the plot is complete. 
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ENDRW must then be called to close the data set. It is 
important to note that in any case, NUM, the data point 
counter must be initialized to 0 at the beginning of each 
run and incremented prior to its respective group of DRWG 
statements . 

It is important "also to remember that ENDRW merely closes out 
the current plot data set and decrements its calling argu- 
ment NPLOT. ENDRW does not initiate a new DSL run. This 
must be done by the user with a call to RERUN, or CONTIN or 
by a parameter study or by any other method given in the 
DSL USER'S MANUAL. Figure A-l is a sample outline of a DSL 
program, showing the use of the subroutines DRWG and ENDRW 
to produce plot output. In Figure A-l, two data sets are to 
be formed. Each data set will contain two plots. The first 
plot in each data set will contain two curves. 

Another important point is that the subroutines DRWG and 
ENDRW are completely independent of the DSL statements 
PREPAR, GRAPH and PRTPLOT. The original IBM system output 
its data to the XY plotter via the statements GRAPH and 
PREPAR. At the Naval Postgraduate School, however, PREPAR, 

GRAPH and PRTPLOT (which at this installation must accompany 
GRAPH) are only used to produce the high-speed printer plots. 

The plot package will handle up to ten plots per run and 
ten curves per plot. Any given curve can have up to 2000 
points. If more than 2000 points are collected for a curve, 
the DSL run will be terminated, and the TERMINAL Region entered. 
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INTGER NUM,NPLOT 
CONST NPLOT=2 

CONTRL FINTIM=10 . 0 , DELT=0 . 1 , DELS=0 . 1 
INITIAL 



NUM=0 

DERIVATIVE 



SAMPLE 

NUM=NUM+1 

CALL DR¥G( 1 , 1 ,NUM, TIME , X) 
CALL DRWG( 1,1 ,NUM, TIME, Y ) 
CALL DRWG(2,1,NUM,X,Y) 



TERMINAL 

CALL ENDRW(NPLOT) 

(parameter or constant changes 
using either FORTRAN assignment 
statements or a DSL END statement 
' followed by the appropriate DSL 
non-executable statements) 

CALL RERUN (when FORTRAN assignments 
have been used) 

END 

STOP 



Figure A-l 



Upon successful completion of the DSL job, the plot 
program PLOTHALL is executed. Each plot set up in the DSL 
job must be accompanied by two title cards. Columns 1 
through ^8 on both cards are reserved for titles. Columns 
75 through 80 of the second title card constitutes the mode 
field (Figure A-2) . 

The mode numbers are listed in Table A-l. If a printed 
output is specified by supplying a printed output interval, 
then a listing of the points actually plotted will be provided 
sampled at the interval specified. The printed output 
interval can have a range from every point plotted to one of 
every 999 points plotted (assuming 999 or more points exist) . 
If column 75 -is left blank, no grid of scaled axes will be 
drawn; just a left-hand and bottom scale will be provided. 

If a 1 is placed in this column, a labeled grid at 3 inch 
intervals will be drawn. This grid is time consuming and 
should be omitted if not essential. 

If a mode other than mode 0 is specified a third title 
card must be provided for the particular plot in question. 

The appropriate data must be entered via this title card 
as shown in Figure A-3- 

Figure A— 4 provides a listing of the JCL (job control 
language) cards necessary to implement a DSL job. 
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J 0 0 0 0 0 0 0 0 0 0 0 C 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 0 : c c c r C no 0 0 r o o 0 0 0 0 0 0 0 0 0 0 0 0 0 C D 0 C 0 0 0 0 0 0 0 0 C 0 0 0 Q 

i 3 < 5 s i > * m? "> '■> ') ‘1~: :• v: :■» : ,f ?t :s t c y -j .• r ^ • •.) : [\ u 5: s ai :o v a 53 to ei 62 a r< «> r v f* 5 * ?j ;t 12 j: :t :: «r. p % ^ 



I 1 M 1 1 1 1 1 1 I 1 1 I 1 1 111 1 1 1 1 1 ! Ill 1 i I! : l 111 1 1i 1 M 1 1 

\ * / , / 

> 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 . 2 2 2 2 7 2 2 2 2 2 2 2 2 2 2 ? 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 ? 2 2 2 2 I 



•444 4 4 4 4444-1 



1 1 U 11 1 1 1 1 1 I 1 1 1 1 1 1 1 1 1 1 11 1 1 1 1 1 l ! 1 

2 2 2 



* 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 .5 3 3 3 3 3 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 

C - •' /V 

4 4 4 4 4 4 4 4 4 4 < - 1 4 - ' ' < 4 ' 4 i i 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 

’ < “> 5 5 5 5 5 5 5 5 5 5 5 5 5 5 r j 5 5 5 5 S 

‘ ' c e f. 6 8 e *■ 



Columns 1 thru 48 on the 
first two title cards per 
plot is the title field. 

The user's name or box 
number should appear 
somewhere within this field 
on either card. 



v 

Mode Field: 

Column 75: grid desired 
Column 76-78: printed output 
interval 

Column 79-80: mode number 



Figure A-2 
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Scaling Fields 



Graph Dimension Fields 



Note: Columns 1 thru ^0 are reserved for Xmin, Dx, 

Ymin and Dy (fields left justified as needed). Values 
may be input in F format or E format. If E format is 
used, data must be right justified in the appropriate 
field. 



Note: Xmin is the value appearing on the extreme left 

end of the horizontal scale and Ymin is the value 
appearing on the bottom of the vertical scale. Dx and 
Dy are the respective scale increments. Lenx and Deny 
are the dimensions of the plot in inches in the X and 
Y directions, where: 



0 < Lenx <_ 9*0 
. 0 < Leny < 20.0 



Obviously Dx and Dy must be compatible with Lenx and 
Leny. When no lengths are specified a 9.0" by 9.0" 
plot will be drawn. 



Note: Columns *J1 thru ^5 and 51 thru 55 are reserved 

for Lenx and Leny (fields left justified as needed). 
Values must be input in F format. 



Figure A- 3 



C 0 0 

v. *J 63 

1 l 1 

2 2 2 
3 3 3 
44 4 
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// (standard green job card) 

// EXEC DSL 
//DSL. INPUT DD * 

(dsl deck goes here) 

//PLOT.SYSIN DD ' x ' (if plots are to be drawn) 
(2 or 3 title cards per plot requested) 



/* 

/« 



(Note: 2 /# ' s must be provided) 



Figure A-4 
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